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In the present work, we examine the combined effects of cubic and quintic terms of the long 
range type in the dynamics of a double well potential. Employing a two-mode approximation, we 
systematically develop two cubic-quintic ordinary differential equations and assess the contributions 
of the long-range interactions in each of the relevant prefactors, gauging how to simplify the ensuing 
dynamical system. Finally, we obtain a reduced canonical description for the conjugate variables of 
relative population imbalance and relative phase between the two wells and proceed to a dynamical 
systems analysis of the resulting pair of ordinary differential equations. While in the case of cubic 
and quintic interactions of the same kind (e.g. both attractive or both repulsive), only a symmetry 
breaking bifurcation can be identified, a remarkable effect that emerges e.g. in the setting of repulsive 
cubic but attractive quintic interactions is a "symmetry restoring" bifurcation. Namely, in addition 
to the supercritical pitchfork that leads to a spontaneous symmetry breaking of the anti-symmetric 
state, there is a subcritical pitchfork that eventually reunites the asymmetric daughter branch with 
the anti-symmetric parent one. The relevant bifurcations, the stability of the branches and their 
dynamical implications are examined both in the reduced (ODE) and in the full (PDE) setting. 



I. INTRODUCTION 

In the study of both atomic and optical physics problems, often analyzed in the realm of nonlinear Schrodinger type 
equations [1, 2 , the study of double well potentials has a prominent position. Such potentials can be straightforwardly 
realized in atomic Bose-Einstein condensates (BECs) through the combination of a parabolic (harmonic) trap with a 
periodic potential. Their experimental realization and subsequent study in BECs with self-repulsive nonlinearity has 
led to numerous interesting observations including tunneling and Josephson oscillations for small numbers of atoms 
in the condensate, and macroscopic quantum self-trapped states for large atom number [3 and symmetry-breaking 
dynamical instabilities |4 . These experimental developments have been accompanied by a larger array of theoretical 
studies on issues such as finite- mode reductions and symmetry-breaking bifurcations [SHE] , quantum effects [13] , and 
nonlinear variants of the potentials [M]. Similar features have also emerged in nonlinear optical settings including 
the formation of asymmetric states in dual-core fibers [15 , self-guided laser beams in Kerr media [16], and optically- 
induced dual-core waveguiding structures in photorefractive crystals [17 . 

On the other hand, a theme that has also been progressively becoming of increasing importance within both of 
these areas of physics is that of long range interactions. In the atomic context, the experimental realization of BECs 
of magnetically polarized ^^Cr atoms [18 (see recent review [19 and for a study of double well effects [20 ), as well 
as the study of dipolar molecules [21], and atoms in which electric moments are induced by a strong external field 
[22^ have been at the center of the effort to appreciate the role of long range effects. On the other hand, in nonlinear 
optics, where nonlocal effects have been argued to be relevant for some time now [23 , numerous striking predictions 
and observations have arisen in the setting of thermal nonlocal media. Among them, we single out the existence of 
stable vortex rings [24 the experimental realization of elliptically shaped spatial solitons [25] and the observation of 
potentially pairwise attracting (instead of repelling as in the standard local cubic media) dark solitons [26] , 

Our aim in the present work is to expand on the framework of studies of double well potentials in the presence of 
nonlocal nonlinear interactions by considering cubic-quintic models. Part of the motivation for doing so consists of the 
fundamental relevance of the cubic-quintic nonlinear Schrodinger. The latter is a model that has been used in a variety 
of physical settings. These include the light propagation in optical media such as non-Kerr crystals [W, chalcogenide 
glasses [28 , organic materials [29 , colloids [30 , dye solutions [31 , and ferroelectrics [32 . It has also been predicted 
that this type of nonlinearity may be synthesized by means of a cascading mechanism [33 . An additional part of 
the motivation stems from an interesting set of observations that were made in an earlier work featuring competing 
cubic nonlinearities, one of which was a cubic local and another was a cubic nonlocal one; see [34 and the discussion 
therein. In that work, it was found that for repulsive nonlocal cubic interactions and attractive local ones, it was 
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possible to tune the prefactors determining the competition so as to produce not only a symmetry breaking, but also 
a symmetry-restoring bifurcation. More recently, a similar conclusion in a local cubic-quintic double well potential 
was reached in [35]. 

Here, we present a framework where the competition of cubic and quintic terms can be systematically quantified. 
In addition, to address the problem from a broader perspective, we consider fully nonlocal interactions both for the 
cubic and the quintic terms, rendering the local case a straightforward special-case scenario of our study. We start our 
presentation of the theoretical analysis of section II by developing a two-mode reduction of the system with both the 
cubic and the quintic terms. We systematically examine all the relevant terms and offer a prescription for assessing 
the dominant contributions to the resulting dynamics of the left and the right well. Following an amplitude-phase 
decomposition and examining the variables associated with the population imbalance of the two wells, and their relative 
phase, we construct the Hamiltonian normal form of the two-mode reduction of the cubic-quintic double well system. 
We then explicitly illustrate how the bifurcation analysis of this normal form encapsulates not only the symmetry 
breaking but also the symmetry restoring. We argue that this cubic-quintic realization is the prototypical one where 
both of these effects can be observed and analytically quantified. Subsequently, in section III, we proceed to test the 
relevant predictions by means of a computational bifurcation analysis, as well as through direct numerical simulations 
(in order to monitor the predicted dynamical instabilities). We find very good agreement with the symmetry breaking 
predictions of the model and even a quite fair agreement with the symmetry restoring ones (which arise in a highly 
nonlinear regime and are hence less amenable to a two-mode analysis). We also quantify the disparity of the analytical 
predictions and numerical results for large values of the nonlocality range parameter. Finally, section IV contains our 
conclusions and some directions for future study. 

II. ANALYTICAL APPROACH FOR THE NLS EQUATION WITH TWO NONLOCAL TERMS 

A. Two- mode approximation 

As indicated above, our fundamental model will be the Id NLS equation in the presence of two nonlocal terms, 
namely the cubic and quintic ones: 

idtilj ^ ^tp = Cilj ^ 8 (^J^ Ri{x - x')\ilj{x')\^dx'^ ^ 6 (^J^ R2{x - x')\tlj{x')\^dx'^ i/j (1) 

with s^S = ±1 and the linear operator will be of the standard Schrodinger type 

C =-{l/2)dl + V{x). 
This encompasses the double-well potential of the form: 

V{x) = {l/2)n^x^ + Vosech^ix/w) 

with Cl being the normalized strength of the parabolic trap and it is <C 1 in a quasi- Id situation in BECs (here the 
effective trap frequency is the ratio of the longitudinal trap strength along the condensate over the one of the tightly 
confined transverse directions). In our study we consider a typical experimentally relevant value of (1 = 0.1, while the 
generally tunable (see e.g. [36 ) parameters of the laser beam forming the light defect are chosen to be Vb = 1 and 
w = 0.5 (which we have found to be fairly typical values representative of the phenomenology to be analyzed below). 
For the kernels R2 we will focus our considerations on either the Gaussian 

1 x^ 
Ri{x) — ^exp(-^) 



or the exponential 



in the spirit of the works of [23]. The key parameter here is the range of the nonlocal interaction parametrized by <j. 
Notice that both kernels in the limit of a ^ tend to a genuinely local interaction (i.e., Ri{x) S{x)). 

We now develop the two-mode approximation in order to obtain a decomposition (or more accurately a Galerkin 
truncation) of the solution over the minimal basis of fundamental states. More specifically, we use an orthonormal 
basis composed by the wave functions {^L,(/>i?} = {('^0 — ui)/^/2^ {uq + i^i)/a/2}, where uq and ui (Fig. 1) are 
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the ground state and the first excited state, respectively, corresponding to the first two eigenvalues of C that are 
cjo = 0.13282 and cji = 0.15571 for our choice of potential parameters above. Notice that these two eigenfunctions 
represent modes with support predominantly on the left and right well, respectively. The eigenfunctions iio,i and the 
rotated basis employed herein of 0L,i? are both shown in Fig. [T] 




FIG. 1: The ground state ixo and first excited state u\ of the potential are shown in the top panel. The rotated orthonormal 
basis of and (\)r (with support, respectively, on the left and right well) is shown in the bottom panels. 



The two-mode approximation is then defined as 

t) = CL{t)(j)L{x) + CR{t)(j)R{x) 



(2) 



where cl and cr are complex time-dependent amplitudes and the approximation consists of the truncation of the 
higher modes within the expansion. Before substituting into the initial GP equation, we notice that the action of the 
linear operator C on our basis elements is as follows: 

C\l) = (ytCL - UJCr)(I)l + {^Cr - UJCL)(t)R 

where Vt = (cjq + cji)/2 and u = {ui — ujq)/2 are linear combinations of the two eigevalues of C respectively to the 
solutions uq^ Ui. Subsequently, substitution of our ansatz of Eq. ([2| in the full nonlinear problem of Eq. ([T]) yields: 

^s\cL\^{cL(t)L ^ CR(j)R) i Ri{x - x')(j)\{x')dx' ^ s\cR\^{cL(t)L ^ CR(j)R) i R2{x - x')(})\{x')dx' 



+5[(cic|^ + |cLpCi^)(^L + (ci4 + CL|ci^p)(^i^] j R{x - x')(^L{x')(^R{x')dx' 

^S\cl\\clH^CrC^r) j R2{x-X')c^\{x')dx' ^8\cr\\clH^CrC^r) j R2{X - X')^%{x')dx' 
-^S [{McL\'\cRfcL + cIc^r' + cI\cl\'cI)^L + {MclIWcr + + C^r\cr\' cD^R • 

■ J R2{x - x')(})l{x')(j?R{x')dx' 

+2(5 [(|cL|'ci4 + \cl\''cr)H + {\cl?\cr\^cl + \cl?cIcI)^r] j R2{x - x')^l{x')^R{x')dx' 
+2(5 [(ci|c^p4 + \clWr\^cr)H + (\cr\^cl + cI\cr\^cI)c^r] j R2{x - x')c^l{x')c^L{x')dx' . 
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In order to project the above equation onto the states <j)L,R we multiply with the respective function (notice that 
the eigenfunctions are real due to the Hermitian nature of the operator C) and integrate. This involves the following 
integrals which will play a fundamental role in our considerations below: 



770 = / / Ri{x — x')(j)\{x')(})\{x)dx' dx^ 

r]2 = J J Ri{x - x')(t)\{x')(l)L{x)(t)R{x)dx'dx, 
7^3 = y j Ri{x - x')(j)L{x')(j)R{x')(l)L{x)(t)R{x)dx'dx, 

from the first nonlocal term, as well as 

r]4 = J J R2{x - x')(j)\{x')(l)\{x)dx' dx, r^g = J J R2{x - x')<j)\{x')<j)\{x')(j)L{x)(j)R{x)dx' dx, 

T]^ = J j R2{x — x')(j)\{x')(j)\{x)dx' dx^ r]g = J j R2{x — x')(p\{x')<pR{x')<p\{x)dx' dx^ 

776 = / / R2{x — x')(j)\{x')(j)L{x)(j)R{x)dx' dx^ rjio = J J R2{x — x')<p\{x')(pR{x')<p\{x)dx' dx^ 

V7 = f f R2{x - x')(l)\{x')(l)\{x')(l)\{x)dx'dx, 7/11 = / / R2{x - x')(t)l{x')(t)R{x')(t)L{x)(l)R{x)dx'dx 



from the second nonlocal term. Some alternatives that are derived if we interchange the variables x and x^ or swap 
L and R can also be equivalently considered. A numerical study of the first four intergrals was already conducted 
in [34], where it was found that typically the integrals 7^2,3 can be considered as negligible in comparison to 7^0 which 
is the dominant term. On the other hand, rji is close to 772,3 for near- local interactions (i.e., for small values of a), 
but becomes comparable to 7^0 as the latter descreases for wide nonlocal interaction ranges (i.e., for large a). The 
criterion that we use to determine whether r]i is negligible or not was rjrei ^ 0.01 where rjrei = ^1 — max(|7/2|, l^sD- 
This yields that r]i remains significant until (i.e., down to) a critical value = 2.96 and 1.56 for the Gaussian and 
exponential kernels, respectively. The dependence of the relevant overlap integrals on the range of the interaction a 
is shown in Fig. |2j 




FIG. 2: The overlap integrals ryo, r/i, 7/2, 773 and 7/4 are shown as a function of the interaction range a for the Gaussian (left) 
and exponential (right) kernels. 

Taking into regard the second nonlocal term (which for simplicity we have assumed to share the same range 
parameter as the first), we can see from Fig. [s] that the integrals 7^5,6,... are always negligible but 7/4 appears to be a 
nontrivial competing term. This is to a certain degree intuitively anticipated, as this represents the dominant term 
associated with the quintic interaction. Adapting the same criterion as in [34 (namely r]rei = ^4 — max(|7^2|, 1^3 1)) 5 we 
incorporate the relevant 7/4 for <j < <Jc = 9.15, 7.01 for the Gaussian and exponential kernel, respectively. According 
to this we may distinguish three cases: 



FIG. 3: The overlap integrals 774, 5,..., n are given here as a function of the interaction range a, for the two kernels in order to 
appreciate the dominance of 774 with respect to the remaining terms for the range a < ctc, where the term with prefactor 774 is 
not negligible with respect to the overall dominant term ryo- 



• The terms 770 and 774 are considered for a < = 2.96 (for the Gaussian kernel); 

• The term with prefactor 771 is added when (T5 < a < (Tc. 

• For cr > (Jc, 774 is omitted and only 770, 771 are taken into account. 

For the first case, the projection of the equation onto the states (I)l,r yields 

icL = {ft - p)cl - ojcr + sr]o\cL\'^CL + Sr]4\cL\'^CL 
ICR = {Q.- fi)cR - ucl + sr]o\cRfcR + Stj^IcrI^cr, 

and by introducing Madelung representation of action-angle or amplitude-phase decomposition {cl^r = PL,i?e*^^'^), 
we obtain 



PL = ^PR sm i 



p - n ^ uj^ cosO - sr]opl -Sr]4pl, 



(3) 



where we have defined the relative phase 6 = Ol — Or and the respective equations for pR and Or can be ob- 
tained by exchanging L and R and using —0 instead of 0. Focusing now on the steady solutions (satisfying 
pL,R = Ol,r = 0), we need to enforce 6> = or tt for non-zero amplitudes. This leads us to symmetric and 

antisymmetric (equal or opposite amplitudes) pairs of solutions, namely for 6> = we have the symmetric (only 

2 

positive ones among the) solutions p\ r = (—sr]o ± V^o ~ ^^V4^{^o — p)) I'^^Va with /i < cjq + 7^, for (5 = — 1 ( 

' V / 4774 

p > ujQ — for S = 1). Also, for ^ = tt, we have (only the positive amplitude ones among) the antisymmetric 

2 2 
solutions p\ R = \ —sr]o ± y^77q — 4^774(6^1 — p) ) /2Sr]4 with p < uji -\- — for S = —1 (resp. p > uji — — ^ for 6 = 1). 

' V / 4?74 47^4 

For the asymmetric solutions one has to solve the polynomial 

^2 

^V^Pl,r + sVopi,R + (^ - p)pI,r + , . = 0, 

which can more conveniently be written as a function of the norm of the solutions (representing the atom number in 
BECs and the optical intensity in optics). Thus, introducing N {N = p\ ^ p'^) yields the quartic polynomial 

(5^7?|7V^ + 3sr]lr]oN'' + {36^7]'^ - r]l{p - n))N^ + (5^7?^ - 2sSr]om{p - - S^u^ - r]^{p -n)=0. 
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For the second case (cr > (75) the integrals 7^0 , Vi^ ^4 taken into account and the projection equations onto the 
states ^L,i?, read: 

icL = {ft- fi)cL - ujcr + scLivolcL]'^ + Vi\cr\'^) + ^V^\cl\^cl 
ICR = {Vt- ii)cr - ujcl + scR{r]Q\cR\^ + + ^r]4\cR\^CR. 

Here, the amphtude-phase decomposition yields 

PL = ^PR sin 6 

'l = /i - + cosO - sr]opl - sr]ip\ - 3r]^p\. 

We can, once again, obtain the set of stationary solutions as follows. When ^ = (symmetric case) the so- 
lutions will be (the positive amplitude ones among) p\ r = ^—5(770 + ^1) i \/ (t^o + ^1)^ — 4(^774 (cjq — p)^ /2(^774 for 

p < uoq ^ + ^1) 5 = —1 {p > ujQ — ^ for (5 = 1) and when = tt (antisymmetric case) the solu- 
4774 4r]4 

tions are (the positive amplitude ones among) p^ j^ = (^—s{r]o -\- r]i) ± \/ {r]o + 771)^ — 4Sr]4{uJi — p)^ /2(57^4 and exist 

for LL < uji -\- ^'^^ for 5 = —1 (ll>uji— ^^-^ for (5 = 1). The asymmetric solutions now, directly in norm 

4774 47/4 
expression, will be given by the polynomial 

S^r]lN^ + {3sr]lr] - sr]lr]i)N^ + {SS^r]^ - r]l{p - - 2Smr]r]i)N^ + 
^{s^T]^ - 2sSr]r]4{p - n) - sr]ir]l)N - Srj^u'^ -rf{p-Q)={) 

with T] here standing for A77 = rjo — rji. 

In the third case, when a > the effect of the quint ic terms is deemed to be negligible and the situation reverts 
to the analysis of ^ and is hence omitted here. 



B. The bifurcation analysis 



In order to derive a more convenient form of the system so that we can proceed to the analysis of the spontaneous 
symmetry breaking (SSB) bifurcation, we introduce the population imbalance between the two wells. 



z={Nl- Nr)/N = (|ci|2 - \cr\^)/N, 



(4) 



where Nl^r = \cl,r\'^ — Pl r ^-nd A'' = Nl+Nr. Together with the relative phase between the two wells 9 = 9l — Or, 
this forms a set of conjugate variables, in which we obtain the dynamical system : 



2uj\/l - z- 
2ujz cos 



This can be written in the Hamiltonian form 



sr]Nz — 5r]iN z. 



_m 

00 

m 

dz 



(5) 



with the Hamiltonian function 



1 



Note that r] stands either for r]o {a < a^) or for Ar] = 770 — ^1 (cr^ < a < dc). The system possesses the stationary 
solutions (critical points) (^1,^1) and (2:2,6^2) with zi = Z2 = 0, Oi = 0^ O2 = tt that correspond to the symmetric 
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and antisymmetric solutions, identified above. Furthermore, the stationary solutions representing the asymmetric 
branches are given by: 



= 0, TT. 



These branches emerge and merge as bifurcations from and to the symmetric or antisymmetric solutions and they 
exist for those values of TV for which > 0. Taking z = 0, we get that 



TV = {-ST] ± ^/n^T8S^)/2Sr]4, N = {- 



-ST] 



8Sr]4Uj)/25r]4. 



(6) 



By substituting {s^5) = (1,-1) or (—1,1) we get the same four possible expressions for A/" as a function of a that 
are displayed in fig. 4 and we denote them with A/'q^, A/'f^, A/^^^ and A^g^ (the subscripts and 2 correspond to the (-) 
signs in the left and right expressions of Eq. (|6|, respectively, while the subscripts 1 and 3 to the (+) signs). One can 
then see that when {s^S) = (1, —1) and demanding that z'^ > 0, one gets that N should either lie in the area outside 



the curves Nq^ and N^^ or in the area inside the curves A/"!^ and N^^. In the case of 



-1 and S = 1, the role of 



the symmetric and anti-symmetric branches gets exchanged in as far as the bifurcation of the asymmetric branch is 
concerned (see also below). 





FIG. 4: The critical values Ng"^ , A^r (left panel) and A^l'', N^"^ (right panel) whenever {s,S) = (1,-1) show when the 
bifurcations appear. More specificaly, the left panel corresponds to the bifurcations that occur on the symmetric branch and 
the right panel for those that occur on the antisymmetric one. 

Importantly, it can be observed in Fig. [4] that A^q^ is always negative, hence it is omitted for the principal case 
considered herein, namely 5 = 1 and S = —1. On the one hand, the critical conclusion of our analysis is that for 
a < 7.52, the system is predicted to have for the anti-symmetric branch both a symmetry breaking bifurcation (at 
N = N2') and a symmetry restoring one that eliminates the asymmetric branch (at N = N^^). On the other hand, 
the right panel suggests that A/"!^, coincide a > 7.52, beyond which there is only a single (symmetry breaking) 
bifurcation. However, as will be discussed below, for large interaction range a this prediction seems to have some 
discrepancy from what actually happens as we will see that in fact, we observe a symmetry restoring bifurcation while 
we do not observe a bifurcation at all in the symmetric branch. To the best of our knowledge, this is the first example 
of an analytical prediction of the existence of a symmetry restoring bifurcation, a feature that is unique to the analysis 
of the normal form of the bifurcation for the cubic-quintic case (and cannot be predicted e.g. in the purely cubic case 
two-mode analysis of [34 ). The new critical points appear or disappear as a pitchfork bifurcation that emerges from 
the antisymmetric solutions for = tt respectively. From the symmetric solution, in this case of 5 = 1 and (5 = — 1, 
only a single bifurcation arises at A^ = A/'f^. 

For the opposite case (to the one principally considered herein) of 5 = — 1 and (5 = 1, i.e., for a focusing cubic 
nonlinearity, the bifurcations emerge from the symmetric branch, while for s = 1, i.e., for a defocusing cubic term, 
then the relevant symmetry breakings arose from the anti- symmetric branch. Thus, in this case, we expect an 
asymmetric branch to bifurcate and break the symmetry at N = A^l^, while it returns to the parent symmetric 
branch restoring the symmetry at A^ = N^^ . On the other hand, for the anti-symmetric waveform with a focusing 
cubic nonlinearity, only a single bifurcation arises at A/" = N^^ . We provide further details of each of these bifurcations 
and their comparison with the full numerics of the underlying NLS model in the next section. 
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From the system of Eqs. ([s]), one can reduce the dynamical evolution to a single second-order ODE: 



z = -Auj'^z - {srjNz + 5r]4^N'^z)\/Aijj'^ - Alo'^z'^ - z^ 
which can also be rewritten in the "position-momentum" variables as: 



\ P = -4w2^ - {sr]Nz (57?4A/'2z) v/4w2 _ 4^2^2 _ j 



(7) 



This renders the system amenable to the phase plane representation of the form shown in Fig. [5j Here we observe 
that there is a stationary solution z = p = which is a fixed point of the center type. However, for the cases when 
{s,5) = (1,-1), for crossing the critical point A^f^ in the case of the symmetric branch and for N G [A/"!^, N^^] in the 

case of the anti-symmetric branch, there appear two more critical points at p = and z = ±W 1 



{srjN ^6r]4N^)^' 

representing the asymmetric solutions. The point (0, 0) is a fixed point of center type before the bifurcation occurs, 
but past the relevant critical number of atoms (or optical intensity), it becomes a saddle as the two new (asymmetric) 
fixed points that appear are of center type. Fig. [5] shows the phase space of the full system, as well as the vicinity of 
the critical points for the Gaussian kernel with cr = 1, N^^ = 4.9862 and N = 5. 




FIG. 5: Top panels: The phase space diagrams of the Hamiltonian system when s = 1 and ^ = — 1 with the Gaussian kernel, 
for cr = 1, A/" = 5, and with Ni^ = 4.9862 (after the new fixed points are created at (dzO.4318, 0)). The left panel displays the 
region of phase space near the symmetric solution (0, 0) (saddle) and the right panel the one near one of the asymmetric fixed 
points (0.4318, 0) (center). The bottom panel shows the full phase space diagram of the system for A" = 5. 
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III. NUMERICAL APPROACH 



A. Stationary solutions 

We now turn to the examination of our analysis against the results of numerical bifurcation analysis (and in the 
next subsection also compare them to direct numerical simulations). We focus here on the case where 5 = 1, J = — 1, 
as we are especially interested in the case of competing interactions; we will briefly also touch upon the case of 5 = — 1 
and ^ = 1. In our numerical computations, the stationary solutions are obtained by using a fixed-point Newton- 
Raphson iteration for a finite difference decomposition of the relevant boundary value problem, with a choice of the 
grid spacing of Ax = 0.1 and employing a parametric (and wherever needed a pseudo-arclength) continuation of the 
solutions with respect to the chemical potential parameter /i (in optics this is the so-called propagation constant). 
The linear stability is analyzed by considering the standard linearization around the stationary solutions tpo in the 
form 

V^(x, t) = V^o + e{a{x)e^^ + 6*(x)e^**). 

This yields the eigenvalue problem 

where the operators are defined as 
1 



/+00 r-\-oo 
K{x - x')\ilJo{x')\^dx' + J / K{x - x')\Mx')fdx' 
-oo ^-oo 
+00 r-\-oo 

+s / K{x - x')tl;o{x)tljQ{x')(l){x')dx' ^26 / K{x - x')tl;o{x)tlJo{x')tljf{x')(l){x')dx' 



and 



L 



/+00 p-\-oo 
K{x-x')Mx')Mx)Hx')dx + 25 / K{x - x')Mx)i^oix')^l{x')(j){x')dx' 
-OO J —OO 



for any real function (j). Instability is guaranteed by the existence of any eigenvalues A of the linearized operator with 
5R(A) 7^ in the sense that perturbations along the corresponding eigendirections will deviate exponentially from the 
corresponding fixed point. Recall that this is also the case for all eigenvalues of our Hamiltonian system, since when 
A is an eigenvalue, so are — A, A"^ and — A"^. In the case where all eigenvalues are found to be purely imaginary, then 
the solution is found to be marginally stable. 

In our specific case of competing interactions, we comment on the following. The positive value (5 = 1) denotes 
the repulsive behavior of the cubic nonlocal term while the negative one S = —1 leads to attractive behavior of the 
quintic nonlocal nonlinearity. As we examine the bifurcation problem of nonlinear states from the corresponding 
linear eigenstates, we expect that for lower values of N (i.e., weaker nonlinearities), the former repulsive term should 
be dominant, while for larger values of N (i.e., stronger nonlinearities), it is anticipated that the latter attractive 
term will take over. This is accurately reflected in the numerical bifurcation diagrams that we now show in Figs. [6][8j 
for three (distinct by roughly an order of magnitude in each case) values of the range a. The first value of a = 0.1 
in Fig. [g] is supposed to refiect the local case (since the range of interaction is much smaller than any other intrinsic 
length scale in the system). Here the agreement with the two-mode approximation is very good quantitatively for 
low N and very good qualitatively (and even good quantitatively for some features such as chemical potentials of 
critical points) for large N. The quality of these types of agreements is found to be preserved for an intermediate 
interaction range of <j = 1 in Fig. [7| However, when the interaction range becomes sufliciently large that it competes 
(or overcomes) the length scale of the potential wells, then fundamental disparities are expected to be found and that 
is the very conclusion of Fig. |8] for a = 8. 

In the first case where a = 0.1, the symmetric and antisymmetric branches of nonlinear states emanate from 
/i = 0.1328 and 0.1557, as expected, respectively {ujq and cji), both of them being dynamically stable, for sufliciently 
small values of N. The rightward bending of the branches for small N confirms the dominance of the self-repulsive 
part of the (cubic) interactions for small A^, as indicated above. The antisymmetric branch (top right panel of Fig. [g] 
and see also the zoom of the bottom panel of the figure) is destabilized and the theoretically predicted asymmetric 
branch emerges. The numerical value of the chemical potential for the bifurcation point is found to be /i = 0.1686, 
whereas the corresponding analytical one is /i = 0.1679, confirming the quantitative nature of the agreement with the 
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FIG. 6: The stationary solution branches for the case s = 1, ^ = — 1 when the interaction range is a — 0.1 expressed in terms 
of the normahzed A/" as a function of ji. The analytical predictions are denoted with the purple dash-dotted line while the 
numerically determined solutions are denoted with the solid line that is blue when it is stable and red otherwise. The top 
left panel shows the symmetric solutions, while the top right presents the antisymmetric ones, both including the asymmetric 
bifurcations that emerge from them. The bottom panel presents a detail of the symmetry-breaking effect, showcasing the 
quality of its approximation by the two-mode expansion. 

two-mode approximation. For larger N ^ we observe that the asymmetric solution has two apparent turning points 
(where the sign of dN/dja changes, but in fact its stability does not change), before it reaches the anti-symmetric 
branch at the numerically computed value /i = 0.381 where we observe the symmetry restoring effect, which, in fact, 
re-stabilizes the anti-symmetric branch. In our theoretical analysis, we observe the same qualitative behavior and the 
symmetry restoring occurs at /i = 0.3723, in reasonable agreement with the full numerical results. Two additional 
observations should be made here. On the one hand, since the symmetry restoring occurs at much larger values of 
the relevant agreement is expected to be less adequate quantitatively than for the symmetry breaking occurring 
at lower N. This is because a two-mode expansion is less appropriate of a reduction at such higher nonlinearities. 
On the other hand, it can indeed be observed that while the overall trend of the two curves is the same (and even 
critical/turning points in terms of their chemical potential are rather accurately captured), this agreement is not 
adequate quantitatively e.g. for critical values of N (or for detailed quantitative matching of the curves for large 
N). For the symmetric solution of the top left panel of Fig. [oj we can observe that it is increasing monotonically 
until /i = 0.359 where it sustains a pitchfork bifurcation leading to the emergence of an asymmetric branch and also 
a subsequent turning point. The symmetric branch becomes unstable thereafter and the asymmetric emerging state 
is the stable daughter branch. Notice that the theoretical analysis is once again quantitatively accurate for small N 
and the agreement becomes more qualitative for higher A^'s. The critical point for the emergence of the asymmetric 
branch is predicted for /i = 0.3492 in reasonable agreement with the full numerical result. 

For the case of a = 1 the effects are similar to those in the previous case. The symmetry breaking of the antisym- 
metric branch (top right, as well as zoom in of the bottom panel of Fig. [7| occurs now at /i = 0.168 according to the 
numerical results and at /i = 0.1673 in the two- mode approximation, again attesting to its validity for small N. After 
following a similar trajectory with the case a = 0.1, the asymmetric solution merges back to the antisymmetric one 
at /i = 0.374 (numerical value) or at /i = 0.364 (analytical value) with the antisymmetric branch again regaining its 
stability past the symmetry restoring bifurcation. The symmetric solution (top left panel of Fig. [tI) again increases 
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FIG. 7: This figure shows the same features as the previous one for the symmetric branch (top left panel), the anti-symmetric 
branch (top right panel) and a zoom-in of the symmetry breaking (bottom panel). However, the interaction range here is an 
order of magnitude larger, namely a = 1. 



monotonically until it sustains a symmetry breaking bifurcation of its own at /i = 0.355. The two- mode approximation 
predicts this bifurcation to arise at /i = 0.342. 

Next, in Fig. |8| we increase the interaction range, roughly, another order of magnitude by setting a = 8. Here, 
as may be intuitively expected given that the interaction range is wider than the wells of the potential, the results 
are quite different. For small values of fi (and thus atom number N or optical power) we have a quite satisfactory 
agreement (even quantititative) with the two mode approximation, as may be expected. As a demonstration of 
that, we note that the symmetry breaking of the antisymmetric branch occurs in our analysis at fi = 0.1981, while 
numerically it is found to take place at /i = 0.195. On the other hand, due to the predicted earlier collision of the 
critical points TV^^ and N^^, there is no symmetry restoring taking place in our normal form reduction. Nevertheless, 
we observe that such a restoring, in fact, still takes place in the full numerical bifurcation diagram. Furthermore, in 
this case, we have not been able to detect a symmetry-breaking bifurcation in the case of the symmetric branch, even 
though such a bifurcation is predicted within the reduction. This illustrates that for such large values of a, even the 
qualitative agreement previously associated with the large N case dynamics should not be expected to be present. 

Finally, we examine also one case where we switch the signs of the nonlocal terms to (s, (^) = (— 1, 1), so now the cubic 
term is the one that behaves attractively while the quintic one behaves repulsively. This is illustrated in Fig. [9j The 
interaction range a is selected here to be 1 and here we see that the same phenomenology appears in a region where the 
cheminal potential varies from —0.08 to 0.155, thus attaining negative values. As earlier, both states emanate for the 
same values of /i and as we decrease its value we observe the symmetry breaking at ja = 0.1212 (both for numerical and 
analytical) this time on the symmetric state which becomes unstable. As we further decrease the chemical potential to 
negative values of /i, the symmetry restoring of the asymmetric state towards its parent symmetric branch occurs at 
/i = —0.0727 (numerical value). The analytical prediction for this critical point is ja = —0.0755. Hence, once again we 
observe a good qualitative agreement for larger N (although once again slight quantitative disparities exist between 
the overall curves and the critical points in terms of N). A look at the antisymmetric branch now shows us that a 
bifurcation occurs at the point where the solution changes slope (dN/d/j.)^ precisely at /i = —0.0465 (numerical) and is 
theoretically predicted to arise at /i = —0.0526 (analytical) with the antisymmetric branch becoming unstable past this 
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FIG. 8: Same as the previous two figures, but now for large nonlocality interaction range in tfie case of cr = 8. 



critical point. Once again the zoom of the bottom panel confirms the quantitative nature of the analytical- numerical 
agreement for small values of which retains its qualitative value even for larger N. 



B. Dynamics 

Finally, we briefly turn to the dynamics of the system, in order to observe the implications of the dynamical 
instability due to the symmetry breaking. The relevant evolution of the unstable solutions for /i = 0.19 and /i = 0.25, 



in the case of a = 1 (recall that 5 = 1 and S = —1) are shown in Fig. 10 In both cases, it can be seen that the weak 
perturbation added on top of the exact numerical solution in the initial conditions has a projection along the unstable 
eigenmode. This projection, for sufficiently long times (about 200 in the left panel and about 100 in the right panel), 
gets amplified and eventually leads to a visible (i.e., of order unity) symmetry breaking in the profile of the state. 
While the space-time evolution of the density (in the atomic case; optical intensity in the optical case) is shown in 



Fig. 10, an interesting alternative way to visualize the instability was proposed recently by ^7\. In the latter work, the 



PDE dynamics was, in fact, projected to the phase plane of the two-mode approximation and visualized therein. An 



example of such a visualization for the case of /i = 0.19 can be seen in Fig. 11 From both the phase plane curves and 
the profiles illustrated underneath of the solution at different times, we can extract some interesting conclusions. In 
particular, in the one degree of freedom reduction of our theoretical analysis, the trajectory occurs over iso-contours 
of the energy. Hence, the kind of phase plane picture shown in Fig. pT] would only be possible by "conglomerating" 
many distinct orbits. However, it is important to appreciate that the PDE has infinitely many degrees of freedom. 
In that capacity, it is possible for the "subspace" of our two-mode approximation to dissipate energy towards (or 
possibly regain energy from) higher energy states (of the point spectrum of the system). In so doing, it appears as if 
the system visits further and further inward trajectories of lower energy, because indeed the excess energy has been 
imparted to other degrees of freedom. This yields a clear illustration of how the subspace of our two-modes is a closed 
system for the ODE reduction, but instead is an open system for the full PDE evolutionary dynamics [41]. 
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FIG. 9: Same as the previous three figures, but now for the focusing cubic/ defocusing quintic case of s = — 1 and ^ = 1, for an 
intermediate interaction range of cr = 1. 



IV. CONCLUSIONS 



In the present work, we examined double well potentials in the presence of nonlocal interactions both in the 
cubic and in the quintic part of the nonlinearity. We attempted to address such settings by means of a two-mode 
decomposition that has the notable advantage that nonlocality is not substantially different to handle therein, as 
the nonlocal kernels merely contribute to relevant overlap integrals that need some systematic book-keeping, but are 
otherwise not considerably harder than is the locally nonlinear case. There are some particularly important attributes 
of the quintic case that we were able to extract via a normal form reduction and phase plane visualization (under 
suitable circumstances of "competition" e.g. for a defocusing cubic but focusing quintic nonlinearity). One such is 
that contrary to the purely cubic case, the reduction is able to predict not only a symmetry breaking bifurcation, 
but also a symmetry restoring one (at least for a suitable interval of range parameters for the interaction kernel). 
Another unusual characteristic is that symmetry breaking bifurcations are encountered both for the symmetric and 
the antisymmetric branch, again differently than is the case for the cubic nonlinearity in the double well setting. 
These features were tested against numerical bifurcation results and good agreement was found where appropriate 
(e.g. low atom numbers and a suitable range of the interaction range). Disparities arising for high N and large a were 
systematically explained. Finally, the instability dynamics was visualized not only by space-time density evolution 
plots but also by offering its projection to the phase plane of the double well theoretical reduction and assessing the 
similarities and differences therein of the ODE approximation and full PDE result. 

There are numerous possibilities for the extension of the present results to more elaborate contexts. On the one 
hand, even in the one-dimensional setting, one could envision a study of different interaction ranges between the 
cubic and quintic terms (or, for that matter, combinations of local and nonlocal nonlinearities within the cubic 
and/or quintic terms). On the other hand, extensions to one dimensional settings with more wells would bring along 
a richer phenomenology (in that setting the three-well local case has been studied [38] and was recently revisited 
in [39 ), while in higher dimensional settings such as 2d, four well settings in a square configuration [40 or other 
configurations exploiting the geometry of the system would be interesting to study. 
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FIG. 10: Spatio-temporal contour plot of the density of the unstable solutions when cr = 1, for s = 1 and ^ = — 1. The panels 
are initialized with (a weakly perturbed case example of) the antisymmetric solution for /i = 0.19 and 0.25 (left and right, 
respectively) . 
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